%
%  Pierre Tremblay, February 2004
%

\defmodule {VarianceGammaProcessDiff}

This class represents a \emph{variance gamma} (VG) process
$\{S(t) = X(t; \theta, \sigma, \nu) : t \geq 0\}$.
This process is generated using \emph{difference of gamma sampling}
(see \cite{fAVR03a,fAVR06a}), which uses the representation of the VG process
as the difference of two independent \class{GammaProcess}'es (see \cite{fMAD98a}):
\begin{equation}
X(t; \theta, \sigma, \nu) := X(0) + \Gamma^{+}(t; \mu_{p}, \nu_{p})
                                  - \Gamma^{-}(t; \mu_{n}, \nu_{n})
\label{dblGammaEqn}
\end{equation}
where $X(0)$ is a constant corresponding to the initial value of the process
and
\begin{equation}
\begin{array}{rcl}
\mu_{p} & = & (\sqrt{ \theta^{2} + 2\sigma^{2}/\nu } + \theta)/2  \\[6pt]
\mu_{n} & = & (\sqrt{ \theta^{2} + 2\sigma^{2}/\nu } - \theta)/2 \\[6pt]
\nu_{p} & = & \nu \mu_{p}^{2}\\[8pt]
\nu_{n} & = & \nu \mu_{n}^{2}
\end{array}
\label{dblGammaParams}
\end{equation}

\bigskip\hrule\bigskip

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{code}
\begin{hide}
/*
 * Class:        VarianceGammaProcessDiff
 * Description:  
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       
 * @since        2004

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */
\end{hide}
package umontreal.iro.lecuyer.stochprocess;\begin{hide}
import umontreal.iro.lecuyer.rng.*;
import umontreal.iro.lecuyer.probdist.*;
import umontreal.iro.lecuyer.randvar.*;

\end{hide}

public class VarianceGammaProcessDiff extends VarianceGammaProcess \begin{hide} {
    protected GammaProcess gpos;
    protected GammaProcess gneg;
    protected double       mup, mun,
                           nup, nun;
\end{hide}
\end{code}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection* {Constructors}
\begin{code}

   public VarianceGammaProcessDiff (double s0, double theta, double sigma,
                                    double nu, RandomStream stream) \begin{hide} {
        this (s0, theta, sigma, nu,
              new GammaProcess (0.0, 1.0, 1.0, stream),
              new GammaProcess (0.0, 1.0, 1.0, stream));
        // Params mu, nu of the 2 gamma processes are redefined in init()
        // which will be called after a call to 'setObservTimes'
    }\end{hide}
\end{code}
\begin{tabb}
Constructs a new \texttt{VarianceGammaProcessDiff} with parameters
$\theta = \texttt{theta}$, $\sigma = \texttt{sigma}$, $\nu = \texttt{nu}$
and initial value $S(t_{0}) = \texttt{s0}$. \texttt{stream} is
used by two instances of \class{GammaProcess}, $\Gamma^{+}$ and $\Gamma^{-}$,
respectively. The other parameters are
as in the class \class{VarianceGammaProcess}.
The \class{GammaProcess} objects for $\Gamma^{+}$ and $\Gamma^{-}$ are
constructed using the parameters from (\ref{dblGammaParams}) and their
initial values  $\Gamma^{+}(t_{0})$ and $\Gamma^{-}(t_{0})$ are set to $0$.
\end{tabb}
\begin{code}

   public VarianceGammaProcessDiff (double s0, double theta, double sigma,
                                    double nu, GammaProcess gpos,
                                    GammaProcess gneg) \begin{hide} {
        this.gpos = gpos;
        this.gneg = gneg;
        setParams (s0, theta, sigma, nu);
        gneg.setStream(gpos.getStream());  // to avoid confusion with stream because there is only
                                 // one stream in the other constructor
    }\end{hide}
\end{code}
\begin{tabb}
The parameters of the \class{GammaProcess}
objects for $\Gamma^{+}$ and $\Gamma^{-}$ are
set to those of (\ref{dblGammaParams}) and their
initial values  $\Gamma^{+}(t_{0})$ and $\Gamma^{-}(t_{0})$ are set to $t_0$.
The \texttt{RandomStream}  of the $\Gamma^{-}$ process is overwritten with the
\texttt{RandomStream}  of the $\Gamma^{+}$ process.
\end{tabb}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection* {Methods}
\begin{code}\begin{hide}

   public double nextObservation() {
        // This implementation takes possible bridge sampling into account
        double s = x0 + gpos.nextObservation() - gneg.nextObservation();
        observationIndex = gpos.getCurrentObservationIndex();
        path[observationIndex] = s;
        observationCounter++;
        return s;
     }

// no longer useful, this method was created to automaticaly alternate
// between the two processes the uniform random variables used in the
// in the simulation.  However, this method does not work if the two
// GammaProcess are PCA...
//    public double[] generatePath()  {
//         gpos.resetStartProcess();
//         gneg.resetStartProcess();
//         double s;
//         for (int i=1; i<=d; i++) {
//            s = x0 + gpos.nextObservation() - gneg.nextObservation();
//            path[gpos.getCurrentObservationIndex()] = s;
//            // Note: we must get the observCounter from gpos in the case that
//            // the process is generated by a Gamma bridge
//            // The observCounter of gneg should be the same because both
//            // gamma processes should be generated the same way
//         }
//         observationIndex   = d;
//         observationCounter = d;
//         return path;
//     }
\end{hide}

   public double[] generatePath() \begin{hide} {
        double[] pathUP = gpos.generatePath();
        double[] pathDOWN = gneg.generatePath();

        for (int i=0; i<d; i++) {
           path[i+1] = x0 + pathUP[i+1] - pathDOWN[i+1];
        }
        observationIndex   = d;
        observationCounter = d;
        return path;
    }\end{hide}
\end{code}
\begin{tabb} Generates, returns and saves the path.  To do so, the
path of $\Gamma^{+}$ is first generated and then the path
of $\Gamma^{-}$.  This is not the optimal way of proceeding
in order to reduce the variance in QMC simulations; for that,
use \texttt{generatePath(double[] uniform01)} instead.
\end{tabb}
\begin{code}

   public double[] generatePath (double[] uniform01) \begin{hide} {
        int dd = uniform01.length;
        int d = dd / 2;

        if (dd % 2 != 0) {
           throw new IllegalArgumentException (
                     "The Array uniform01 must have a even length");
        }

        double[] QMCpointsUP = new double[d];
        double[] QMCpointsDW = new double[d];

        for(int i = 0; i < d; i++){
            QMCpointsUP[i] = uniform01[2*i];  // keeps the odd numbers for the gamma process
            QMCpointsDW[i] = uniform01[2*i + 1]; // and the even for the BM process
        }
        gpos.resetStartProcess();
        gneg.resetStartProcess();

        double[] pathUP = gpos.generatePath(QMCpointsUP);
        double[] pathDOWN = gneg.generatePath(QMCpointsDW);

        for (int i=0; i<d; i++) {
           path[i+1] = x0 + pathUP[i+1] - pathDOWN[i+1];
        }
        observationIndex   = d;
        observationCounter = d;
        return path;
    }
\end{hide}
\end{code}
\begin{tabb}
Similar to the usual \texttt{generatePath()}, but here the uniform
random numbers used for the simulation must be provided to the method.  This
allows to properly use the uniform random variates in QMC simulations.
This method divides the table of uniform random
numbers \texttt{uniform01} in two smaller tables, the first one containing
the odd indices of \texttt{uniform01} are used to generate the path of $\Gamma^{+}$
and the even indices are used to generate the path of  $\Gamma^{-}$.
This way of proceeding further reduces the
variance for QMC simulations.
\end{tabb}
\begin{code}

   public void resetStartProcess() \begin{hide} {
        observationIndex   = 0;
        observationCounter = 0;
        gpos.resetStartProcess();
        gneg.resetStartProcess();
    }\end{hide}
\end{code}
\begin{tabb} Sets the observation times on the \texttt{VarianceGammaProcessDiff}
as usual, but also applies the \texttt{resetStartProcess} method to the two
\class{GammaProcess} objects used to generate this process.
\end{tabb}
\begin{code}

   public GammaProcess getGpos() \begin{hide} {
        return gpos;
    }\end{hide}
\end{code}
\begin{tabb} Returns a reference to the \class{GammaProcess} object \texttt{gpos}
used to generate the $\Gamma^{+}$ component of the process.
\end{tabb}
\begin{code}

   public GammaProcess getGneg() \begin{hide} {
        return gneg;
    }\end{hide}
\end{code}
\begin{tabb} Returns a reference to the \class{GammaProcess} object \texttt{gneg}
used to generate the $\Gamma^{-}$ component of the process.
\end{tabb}
\begin{code}\begin{hide}

    protected void init() {
        // super.init() is not called because the init() in VarianceGammaProcess
        // is very different.
        mup = 0.5 * (Math.sqrt (theta*theta + 2*sigma*sigma/nu) + theta);
        mun = 0.5 * (Math.sqrt (theta*theta + 2*sigma*sigma/nu) - theta);
        nup = mup * mup * nu;
        nun = mun * mun * nu;
        if (observationTimesSet) {
            path[0] = x0;
            gpos.setParams(t[0], mup, nup);
            gneg.setParams(t[0], mun, nun);
        }
    }\end{hide}

   public void setObservationTimes (double t[], int d) \begin{hide} {
         gpos.setObservationTimes(t, d);
         gneg.setObservationTimes(t, d);
         super.setObservationTimes(t, d);
// the initial value is set to t[0] in the init, which is called in
// super.setObservationTimes(t, d).
     }\end{hide}
\end{code}
\begin{tabb} Sets the observation times on the \texttt{VarianceGammaProcesDiff}
as usual, but also sets the observation times of the underlying
\class{GammaProcess}'es.
\end{tabb}
\begin{code}

   public RandomStream getStream() \begin{hide} {
      return gpos.getStream ();
   }\end{hide}
\end{code}
\begin{tabb}
Returns the \texttt{RandomStream} of  the $\Gamma^{+}$ process.
\end{tabb}
\begin{code}

   public void setStream (RandomStream stream) \begin{hide} {
         gpos.setStream(stream);
         gneg.setStream(stream);
   }\end{hide}
\end{code}
\begin{tabb}
Sets the \externalclass{umontreal.iro.lecuyer.rng}{RandomStream}
of the two \class{GammaProcess}'es to  \texttt{stream}.
\end{tabb}
\begin{code}\begin{hide}
}
\end{hide}\end{code}
